home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / MTDS.FOR < prev    next >
Text File  |  1985-11-29  |  5KB  |  155 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE MTDS
  5. C
  6. C        PURPOSE
  7. C           MULTIPLY A GENERAL MATRIX A ON THE LEFT OR RIGHT BY
  8. C           INVERSE(T),INVERSE(TRANSPOSE(T)) OR INVERSE(TRANSPOSE(T*T))
  9. C           THE TRIANGULAR MATRIX T IS STORED COLUMNWISE IN COMPRESSED
  10. C           FORM, I.E. UPPER TRIANGULAR PART ONLY.
  11. C
  12. C        USAGE
  13. C           CALL MTDS(A,M,N,T,IOP,IER)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           A     - GIVEN GENERAL MATRIX WHITH M ROWS AND N COLUMNS.
  17. C           M     - NUMBER OF ROWS OF MATRIX A
  18. C           N     - NUMBER OF COLUMNS OF MATRIX A
  19. C           T     - GIVEN TRIANGULAR MATRIX STORED COLUMNWISE UPPER
  20. C                   TRIANGULAR PART ONLY. ITS NUMBER OF ROWS AND
  21. C                   COLUMNS K IS IMPLIED BY COMPATIBILITY.
  22. C                   K = M IF IOP IS POSITIVE,
  23. C                   K = N IF IOP IS NEGATIVE.
  24. C                   T OCCUPIES K*(K+1)/2 STORAGE POSITIONS.
  25. C           IOP   - INPUT VARIABLE FOR SELECTION OF OPERATION
  26. C                   IOP = 1 - A IS REPLACED BY INVERSE(T)*A
  27. C                   IOP =-1 - A IS REPLACED BY A*INVERSE(T)
  28. C                   IOP = 2 - A IS REPLACED BY INVERSE(TRANSPOSE(T))*A
  29. C                   IOP =-2 - A IS REPLACED BY A*INVERSE(TRANSPOSE(T))
  30. C                   IOP = 3 - A IS REPLACED BY INVERSE(TRANSPOSE(T)*T)*A
  31. C                   IOP =-3 - A IS REPLACED BY A*INVERSE(TRANSPOSE(T)*T)
  32. C           IER   - RESULTING ERROR PARAMETER
  33. C                   IER =-1 MEANS M AND N ARE NOT BOTH POSITIVE
  34. C                                 AND/OR IOP IS ILLEGAL
  35. C                   IER = 0 MEANS OPERATION WAS SUCCESSFUL
  36. C                   IER = 1 MEANS TRIANGULAR MATRIX T IS SINGULAR
  37. C
  38. C        REMARKS
  39. C           SUBROUTINE MTDS MAY BE USED TO CALCULATE THE SOLUTION OF
  40. C           A SYSTEM OF EQUATIONS WITH SYMMETRIC POSITIVE DEFINITE
  41. C           COEFFICIENT MATRIX. THE FIRST STEP TOWARDS THE SOLUTION
  42. C           IS TRIANGULAR FACTORIZATION BY MEANS OF MFSD, THE SECOND
  43. C           STEP IS APPLICATION OF MTDS.
  44. C           SUBROUTINES MFSD AND MTDS MAY BE USED IN ORDER TO CALCULATE
  45. C           THE PRODUCT TRANSPOSE(A)*INVERSE(B)*A WITH GIVEN SYMMETRIC
  46. C           POSITIVE DEFINITE B AND GIVEN A EFFICIENTLY IN THREE STEPS
  47. C           1) TRIANGULAR FACTORIZATION OF B (B=TRANSPOSE(T)*T)
  48. C           2) MULTIPLICATION OF A ON THE LEFT BY INVERSE(TRANSPOSE(T))
  49. C              A IS REPLACED BY C=INVERSE(TRANSPOSE(T))*A
  50. C           3) CALCULATION OF THE RESULT FORMING TRANSPOSE(C)*C
  51. C
  52. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  53. C           NONE
  54. C
  55. C        METHOD
  56. C           CALCULATION OF X = INVERSE(T)*A IS DONE USING BACKWARD
  57. C           SUBSTITUTION TO OBTAIN X FROM T*X = A.
  58. C           CALCULATION OF Y = INVERSE(TRANSPOSE(T))*A IS DONE USING
  59. C           FORWARD SUBSTITUTION TO OBTAIN Y FROM TRANSPOSE(T)*Y = A.
  60. C           CALCULATION OF Z = INVERSE(TRANSPOSE(T)*T)*A IS DONE
  61. C           SOLVING FIRST TRANSPOSE(T)*Y = A AND THEN T*Z = Y, IE.
  62. C           USING THE ABOVE TWO STEPS IN REVERSE ORDER
  63. C
  64. C     ..................................................................
  65. C
  66.       SUBROUTINE MTDS(A,M,N,T,IOP,IER)
  67. C
  68. C
  69.       DIMENSION A(1),T(1)
  70.       DOUBLE PRECISION DSUM
  71. C
  72. C        TEST OF DIMENSION
  73.       IF(M)2,2,1
  74.     1 IF(N)2,2,4
  75. C
  76. C        ERROR RETURN IN CASE OF ILLEGAL DIMENSIONS
  77.     2 IER=-1
  78.       RETURN
  79. C
  80. C        ERROR RETURN IN CASE OF SINGULAR MATRIX T
  81.     3 IER=1
  82.       RETURN
  83. C
  84. C        INITIALIZE DIVISION PROCESS
  85.     4 MN=M*N
  86.       MM=M*(M+1)/2
  87.       MM1=M-1
  88.       IER=0
  89.       ICS=M
  90.       IRS=1
  91.       IMEND=M
  92. C
  93. C        TEST SPECIFIED OPERATION
  94.       IF(IOP)5,2,6
  95.     5 MM=N*(N+1)/2
  96.       MM1=N-1
  97.       IRS=M
  98.       ICS=1
  99.       IMEND=MN-M+1
  100.       MN=M
  101.     6 IOPE=MOD(IOP+3,3)
  102.       IF(IABS(IOP)-3)7,7,2
  103.     7 IF(IOPE-1)8,18,8
  104. C
  105. C        INITIALIZE SOLUTION OF TRANSPOSE(T)*X = A
  106.     8 MEND=1
  107.       LLD=IRS
  108.       MSTA=1
  109.       MDEL=1
  110.       MX=1
  111.       LD=1
  112.       LX=0
  113. C
  114. C        TEST FOR NONZERO DIAGONAL TERM IN T
  115.     9 IF(T(MSTA))10,3,10
  116.    10 DO 11 I=MEND,MN,ICS
  117.    11 A(I)=A(I)/DBLE(T(MSTA))
  118. C
  119. C        IS M EQUAL 1
  120.       IF(MM1)2,15,12
  121.    12 DO 14 J=1,MM1
  122.       MSTA=MSTA+MDEL
  123.       MDEL=MDEL+MX
  124.       DO 14 I=MEND,MN,ICS
  125.       DSUM=0.D0
  126.       L=MSTA
  127.       LDX=LD
  128.       LL=I
  129.       DO 13 K=1,J
  130.       DSUM=DSUM-T(L)*A(LL)
  131.       LL=LL+LLD
  132.       L=L+LDX
  133.    13 LDX=LDX+LX
  134.       IF(T(L))14,3,14
  135.    14 A(LL)=(DSUM+A(LL))/T(L)
  136. C
  137. C        TEST END OF OPERATION
  138.    15 IF(IER)16,17,16
  139.    16 IER=0
  140.       RETURN
  141.    17 IF(IOPE)18,18,16
  142. C
  143. C        INITIALIZE SOLUTION OF T*X = A
  144.    18 IER=1
  145.       MEND=IMEND
  146.       MN=M*N
  147.       LLD=-IRS
  148.       MSTA=MM
  149.       MDEL=-1
  150.       MX=0
  151.       LD=-MM1
  152.       LX=1
  153.       GOTO 9
  154.       END
  155.